rm(list = ls())
########################################
reference_system<-'+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs '
xmin=-10136161
xmax =-10077458
ymin= 1800245
ymax=1858949

########################################
library(parallel)
library(pacman)
p_load(dplyr,sf,raster,stars,tidyverse)

if(Sys.info()[['user']]=="santiago.fernandez"){
  if(Sys.info()[["sysname"]]=="Windows"){
    main<- 'C:/Users/santi_2wcovuj/Desktop/Predictions/forest_predictions/CreatedData'
  }
  if(Sys.info()[["sysname"]]=="Linux"){
    main<- '/dataeco01/economia/proyectos/santisap/Deforestation/CreatedData'
  }}else if(Sys.info()[["user"]]=="mineria.ra"){
  main<- 'C:/Users/Mineria.ra/Dropbox/Research/PES Chiapas RCT/CreatedData'
  }

#Function to create id by pixel 

coord_to_smartchiapas<-function(x_vec, y_vec){
  
  # X and Y minimum values of Colombia in input coordinate system

  x_min<- -10136161
  y_min<- 1800245
  #side of the grid square (maxx-minx) (540 in CoMiMo)
  grid_size <- 5
  
  ## Trick vector filled with minimum values for x and y
  x_min_vec<-matrix(x_min, nrow=length(x_vec))
  y_min_vec<-matrix(y_min, nrow=length(y_vec))
  
  ## Vector filled with the id_grid
  #Multiply by 10000 so first 4 coordinates are x and second 4 y
  #Calculate distance in N of grids to the left of Colombia
  #as.integer((x_vec-x_min_vec)/grid_size)
  #add 1000 to avoid codes starting with 0 and keep 8 digits
  id_vec<-
    10000*(as.integer((x_vec-x_min_vec)/grid_size)+1000)+
    ((as.integer((y_vec-y_min_vec)/grid_size)+1000)+1000)
  
  return(id_vec)
}

setwd(main)

raster_path<- "pixel_id/Rasters/"
files<-list.files(raster_path,pattern = '.tif')


columns <- c("August2022" ,"July2021" ,"July2022", "June2020","May2021",
             "May2022"   , "Sept2022"  , "Property" , "Contract","Conafor")


create_pixel_shape<-function(path){
  
  raster <- brick(paste0(raster_path
                       ,path))
  

  rast <- stars::st_as_stars(raster)
  

  
  shape<-st_as_sf(rast)%>%
    st_centroid()


  separated_coord <- shape %>%
    mutate(long = unlist(map(shape$geometry,1)),
           lat = unlist(map(shape$geometry,2)))
  
  
  final_shape<-separated_coord[-c(1:5)]%>%
    mutate(
      id_grid = coord_to_smartchiapas(long,lat)
    )%>%
    mutate(
      polygon_id=gsub("pol_","",gsub(".tif","",path))
    )%>%
    
    st_transform(crs="EPSG:4326")

  
  colnames(final_shape)[1:length(columns)] <- columns
  
  
  final_shape <- subset(final_shape, select = -c(long, lat))#%>%

  write_sf(final_shape,
           paste0("pixel_id/ByPixel_DataResamp/",gsub(".tif",".shp",path)),
           overwrite= T)
  csv <- final_shape%>%
    st_drop_geometry()
  
  write.csv(csv,paste0("pixel_id/ByPixel_DataResamp/",gsub(".tif",".csv",path)),
            row.names = F)
}



# to_Ha <- 4.588438*4.560103/10000


numCores <- detectCores()-6
cl <- makeCluster(numCores) 
clusterExport(cl, ls())
clusterEvalQ(cl,library(stars))
clusterEvalQ(cl,library(tidyverse))
clusterEvalQ(cl,library(sf))
clusterEvalQ(cl,library(raster))

parLapply(cl,files,create_pixel_shape)
stopCluster(cl)


